# Calculate fixed effects for each individual
model_fixef <- feols_custom(
  r2_choose_comparator ~ pair_includes_trans_alt:factor(ind_id) + item_diff_value + r2_reliability_diff + r2_reliability_shown + video_type + comparator_order_in_pair + delivery_incentive_exp + phase + factor(round),
  data = r2_choices_with_item_val %>% filter(discuss_type %in% c("control", "discussion_full")),
  cluster = "group_id"
)

# r2_choices_with_item_val$pair
coeff_items <- model_fixef %>% get_coeff("item_diff_value")

fixefs <- model_fixef %>% tidy() %>%
  filter(str_detect(term, "pair_includes_trans_alt\\:factor\\(ind_id\\)")) %>%
  mutate(ind_id = term %>% str_replace("pair_includes_trans_alt\\:factor\\(ind_id\\)", "")) %>%
  mutate(estimate = estimate / coeff_items) %>%
  mutate(group_id = ind_id %>% str_replace("_.*", "")) %>%
  mutate(group_mean = mean_na(estimate), .by = group_id) %>%
  mutate(group_resid = estimate - group_mean) %>%
  left_join(df %>% select(ind_id, discuss_type))

fixefs %>% group_by(discuss_type) %>% summarise(mean_na(estimate))


r2_choices_with_item_val %>% count_prop(item_diff_value)

# Calculate statistics by treatment group
stats_by_group <- fixefs %>%
  group_by(discuss_type) %>%
  summarize(
    mean_estimate = sprintf("%.1f", mean(estimate, na.rm = TRUE)),
    prob_positive = sprintf("%.2f", mean(estimate >= 0, na.rm = TRUE))
  ) %>%
  ungroup()

  
stats_by_group %>% filter(discuss_type=="control") %>% pull(mean_estimate) %>% as.numeric() %>% write_stat("outputs/stats/pref_mean_control.tex")
stats_by_group %>% filter(discuss_type=="discussion_full") %>% pull(mean_estimate) %>% as.numeric() %>% write_stat("outputs/stats/pref_mean_discussion.tex")

stats_by_group %>% filter(discuss_type=="control") %>% pull(prob_positive) %>% as.numeric() %>% write_percentage("outputs/stats/pref_prob_control.tex")
stats_by_group %>% filter(discuss_type=="discussion_full") %>% pull(prob_positive) %>% as.numeric() %>% write_percentage("outputs/stats/pref_prob_discussion.tex")

# Create the plot with annotations
fixefs %>%
  ggplot(aes(x = estimate, fill = discuss_type)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, alpha = 0.8, 
                 color = "white", linewidth = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50", size = 0.5) +
  scale_fill_brewer(palette = "Set2", direction = -1) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
  labs(x = "WTP to select transgender worker\n(Individual-level, Rs.)", y = "Density") +
  facet_wrap(~discuss_type, labeller = labeller(
    discuss_type = c("control" = "No discussion (private)", "discussion_full" = "3-person discussion")
  )) +
  # Add text annotations from the stats dataframe
  geom_text(
    data = stats_by_group,
    aes(
      x = Inf, y = Inf,
      label = paste0("Mean = ", mean_estimate, "\n", "P(WTP ≥ 0) = ", prob_positive)
    ),
    hjust = 1.1, vjust = 1.2,
    size = 3.5
  ) +
  theme_bw() +
  theme(
    legend.position = "none",
    panel.grid = element_blank(),
    strip.background = element_rect(fill = "white", color = "gray80"),
    strip.text = element_text(face = "bold"),
    axis.title = element_text(face = "bold"),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    plot.margin = margin(t = 5, r = 15, b = 15, l = 5, unit = "pt")
  )

ggsave("outputs/figs/pref_distribution.pdf", width = 7, height = 4)